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CRASHa: coupling continuum and line radiative transfer. 
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1 INTRODUCTION 



The detection of Lya lines from local and distant objects has always 
been of great importance in astrophysics. It has been extensively 
used as indicator of redshift, as a measurement of the star formation 
activity of galaxies and as a probe of their internal structure. In 
the last few years an increasing interest has been devoted to the 
search of Lya emitters (LAEs) at high redshift, which are expected 
to be characterized by a strong Lya emission (Partridge and Peebles 
1968), but significantly attenuated by dust absorption. In fact, it 
has been necessary to wait for dedicated large programs of deep 
narrow band searches like the Large Area Lyman Alpha (LALA) 
and the Subaru Deep Field survey to detect a significant number 
of emission galaxies at high redshift (e.g. Stern et al. 2005; lye 
et al. 2006) and to get complete spectroscopic samples of LAEs 
at redshift z = 4.5, z = 5.7 and z = 6.5 (e.g. Hu et al. 1998, 
2004; Rhoads & Malhotra 2001; Kodaira et al. 2003; Taniguchi et 
al. 2005; Kashikawa et al. 2006; Murayama et al. 2007; Dawson 
et al. 2007). Strong lensing magnification has been necessary to 
move the detection frontiers even further, with several candidates 



ABSTRACT 

In this paper we present CRASHa , the first radiative transfer code for cosmological applica- 
tion that follows the parallel propagation of Lya and ionizing photons. CRASHa is a version 
of the continuum radiative transfer code CRASH with a new algorithm to follow the propaga- 
tion of Lya photons through a gas configuration whose ionization structure is evolving. The 
implementation introduces the time evolution for Lya photons (a feature commonly neglected 
in line radiative transfer codes) and, to reduce the computational time needed to follow each 
scattering, adopts a statistical approach to the Lya treatment by making extensive use of pre- 
compiled tables. These tables describe the physical characteristics of a photon escaping from 
a gas cell where it was trapped by scattering as a function of the gas temperature/density and 
of the incoming photon frequency. With this statistical approach we experience a drastic in- 
crease of the computational speed and, at the same time, an excellent agreement with the full 
Lya radiative transfer computations of the code MCLya. We find that the emerging spectra 
keep memory of the ionization history which generates a given ionization configuration of the 
gas and, to properly account for this effect, a self-consistent joint evolution of line and ioniz- 
ing continuum radiation as implemented in CRASHa is necessary. A comparison between the 
results from our code and from Lya scattering alone on a fixed HI density field shows that the 
extent of the difference between the emerging spectra depends on the particular configuration 
considered, but it can be substantial and can thus affect the physical interpretation of the prob- 
lem at hand. These differences should furthermore be taken into account when computing the 
impact of the Lya radiation on e.g. the observability of the 21 cm line from neutral hydrogen 
at epochs preceeding complete reionization. 
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currently observed up to z ~ 10 (e.g. Pello et al. 2007; Stark et al. 
2007). 

The intense activity reported above is explained by the great 
interest in using LAEs as cosmological probes: LAEs are in fact 
the objects with the highest known z and can be used to study 
large scale structures and galaxy formation in the high redshift uni- 
verse. Number counts, together with the statistics of line shapes, 
are extremely powerful observables from which inferring important 
information about e.g. the properties of the intergalactic medium 
(IGM) and the reionization era (e.g. Hu 2002; Rhoads 2003; Stern 
2005), the photoionization processes and UV photon production 
(e.g. Stark 2007), the tomography of neutral gas, gas velocity field 
and star formation activity (e.g. Kodaira et al. 2003). This predict- 
ing power relies on the fact that the Lya line is shaped inside the 
galaxy interstellar medium and at high redshift is also affected by 
the IGM opacity, which becomes non negligible to Lya photons at 
z > 6 (Fan et al. 2006). 

Emission of Lya photons from high redshift sources has also 
an impact on the detectability of 21 cm line from neutral hydro- 
gen in the IGM. At high redshift, the Wouthuysen-Field effect 
(Wouthuysen 1952; Field 1958, 1959) is in fact extremely efficient 
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in decoupling the spin temperature of the gas, T„, from the cosmic 
microwave background (CMB) temperature, Tomb, allowing the 
21 cm signal to be visible either in absorption or in emission. Fluc- 
tuations in the Lya flux, due both to inhomogeneous distribution 
of the Lya radiation sources and to the scattering in the wings, can 
modify the expected signal (e.g. Barkana & Loeb 2005; Chen & 
Miralda-Escude 2006; Chuzhoy & Zheng 2007; Semelin, Combes 
& Baek 2007), at least as long as a strong Lya background is not es- 
tablished and the radiation intensity reaches a saturation level (e.g. 
Ciardi & Madau 2003; Ciardi & Salvaterra 2007). For these rea- 
sons, it is important to follow the propagation of Lya photons rather 
than assume a homogeneous background as it is generally done. 

Due to the resonant nature of the Lya line propagation, a 
self-consistent and detailed treatment of the line radiation trans- 
fer is required in order to model properly how Lya radiation af- 
fects the IGM, as well as to understand how different physical 
processes shape the spectral features of LAEs. As a consequence 
of the great interest in this field, several semi-analytic and nu- 
merical studies of the Lya radiative transfer have followed the 
first pioneering papers on the subject (Osterbrock 1962; Avery 
& House 1968; Adams 1972; Harrington 1973; Neufeld 1990). 
Analytic solutions have been derived only for few simple geo- 
metrical gas configurations: static plane parallel slabs including 
dust (Neufeld 1990), static uniform sphere (Dijkstra, Haiman & 
Spaans 2006) and uniform gas with pure Hubble flow around a 
steady Lya source (Loeb & Rybicki 1999). Given the difficulties 
in the treatment of radiative transfer though, also the numerical 
approaches developed so far, mostly based on Monte Carlo tech- 
niques, have been in most cases specifically designed for particular 
physical configurations and problems: ID dusty and optically thick 
media (Ahn, Lee & Lee 2000, 2001); 3D arbitrary distribution of 
dustless gas with arbitrary bulk velocity field (Zheng & Miralda- 
Escude 2002); spherically symmetric collapsing gas clouds (Dijk- 
stra, Haiman & Spaans 2006); Lya scattering off opaque, dusty and 
moving clouds (Hansen & Oh 2006); Hubble like expansion flows 
of neutral gas (Loeb & Rybichi 1999; Kobayashi & Kamara 2004). 
Other codes have been specifically designed for studying LAEs and 
Lya pumping in a cosmological context (Gould & Weinberg 1996; 
Cantalupo et al. 2005; Tasitsiomi 2006; Semelin, Combes & Baek 
2007). Verhamme, Schaerer & Maselli (2006, VSM06) have devel- 
oped a general-purpose 3D Lya radiation transfer code applicable 
to dusty media with arbitrary geometries and velocity fields. 

So far analytical, semi-analytical as well as numerical stud- 
ies perform the Lya radiative transfer as a post-process calculation 
by assuming a fixed ionization structure of the gas through which 
it propagates, while none of them has tackled the Lya radiative 
transfer problem by taking into account the effect of an evolving 
ionization configuration. Nevertheless, as we show in the follow- 
ing of this paper, this approximation results to be a poor one for 
some applications of interest, in particular for cosmological stud- 
ies at high redshift, but also when modeling the Lya emission from 
young galaxies. 

In this paper we present CRASHa, a new radiative trans- 
fer scheme which, for the first time in the literature, follows si- 
multaneously the propagation of Lya and ionizing radiation self- 
consistently. This allows us to investigate the effects of evolv- 
ing ionization configurations on the propagation of Lya radiation 
and on the shaping of the line emerging from single objects. The 
impact of an evolving ionization structure can in fact be signif- 
icant and needs to be taken into account: the large cross-section 
of Lya photons makes propagation dominated by resonant scatter- 
ing with HI atoms and the random-walk-like nature of the process 



makes the characteristic time for Lya photon propagation much 
larger than the one for ionizing radiation. If an ionizing continuum 
changes the ionization of the gas through which the Lya photon is 
propagating, the amount of scattering suffered by the line photons 
before escaping will depend on the ionization history of the system. 
In this case, a joint treatment of both line and continuum transfer 
is needed to study the alterations in the Lya spectrum occurring 
during the evolutionary stages of the ionized regions. 

The code presented in this paper is the first step in this direc- 
tion. CRASHa has been implemented as an extension of the 3D ray- 
tracing radiative transfer code for ionizing radiation CRASH (Cia- 
rdi et al. 2001; Maselli, Ferrara & Ciardi 2003; Maselli & Ferrara 
2005; Maselli, Ciardi & Kanekar 2008), by developing a new inde- 
pendent algorithm which follows the path of line photons in time 
and space. As described in details in the following, this new al- 
gorithm makes extensive use of pre-compiled tables which have 
been derived by using the line transfer code MCLya (VSM06) and 
allows to compute in an extremely efficient way the path of line 
photons in arbitrary 3D gas distributions. 

The paper is structured as follows. Section 2 is dedicated to a 
brief overview of CRASH and MCLya, while in Section 3 we de- 
scribe the new method. Some validation tests are shown in Section 
4. In the last Section we present a summary of the paper. 



2 CRASH AND MCLYa 

In this Section we give a brief description of the codes CRASH 
and MCLya for the sake of providing a proper background for 
the description of CRASHa given in Section 3. A more detailed de- 
scription of the two codes is already in the literature: the details 
of CRASH implementation are given mostly in Maselli, Ferrara & 
Ciardi (2003), with updates on a new scheme for the background 
radiation field given in Maselli & Ferrara (2005) and on the latest 
version of the code in Maselli, Ciardi & Kanekar (2008). MCLya 
algorithm is fully described in VSM06. Note that some nomencla- 
ture has been changed for clarity. 



2.1 CRASH: continuum radiative transfer 

CRASH is a 3D ray-tracing radiative transfer code based on Monte 
Carlo (MC) techniques that are used to sample the probability dis- 
tribution functions (PDFs) of several quantities involved in the cal- 
culation, e.g. spectrum of the sources, emission direction, optical 
depth. The MC approach and the code architecture assure a great 
flexibility in the application to a wide range of astrophysical prob- 
lems and allow additional physics to be easily added with a mini- 
mum effort. 

The algorithm follows the propagation of the ionizing radia- 
tion through an arbitrary H/He static density field and at the same 
time computes the variations in temperature and ionization state of 
the gas. Both multiple point sources, located arbitrarily in the box, 
and diffuse radiation (e.g. the ultraviolet background or the radia- 
tion produced by H/He recombinations) can be accounted for. In 
this paper we neglect the treatment of any background radiation for 
simplicity. 

The energy emitted by point sources in ionizing radiation is 
discretized into photon packets, beams of ionizing photons, emit- 
ted at regularly spaced time intervals. More specifically, the total 
energy radiated by a single source of luminosity L s , during the total 
simulation time, t sim , is E s = j " m L s (t s )dt s . For each source, 
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E a is distributed in N p photon packets, emitted at the source lo- 
cation at regularly spaced time intervals, dt — t s im/N p . The time 
resolution of a given run is thus fixed by N p and the time evolution 
is marked by the packets emission: the j-th packet is emitted at time 
tim,e = J x dt, with j — 0, (N p — 1). Thus, the total number 



of emissions of continuum photon packets is N e - 



N v . In its 



latest version (Maselli, Ciardi & Kanekar 2008), the code allows 
for polychromatic packets whose content consists of photons dis- 
tributed in various frequency bins which are populated according 
to the spectral shape assigned to the source. 

The emission direction of each photon packet is assigned by 
MC sampling the angular PDF characteristic of the source. The 
propagation of the packet through the given density field is then 
followed and the impact of radiation-matter interaction on the gas 
properties is computed on the fly. Each time the packet pierces a 
cell i, the cell optical depth for ionizing continuum radiation, rl, 
is estimated summing up the contribution of the different absorbers 
(HI, Hel, Hell). As the probability for a single photon to be ab- 
sorbed in such a cell is: 



P(rt) = l-e 



(1) 



the number of photons absorbed in the cell i is the fraction P(t£) 
of packet content when entering the cell. In the polychromatic im- 
plementation, the same argument applies to the number of photons 
contained in each single frequency bin. The trajectory of the packet 
is followed until its photon content is extinguished or, if continuum 
boundary conditions are not assumed, until it exits the simulation 
volume. 

The time evolution of the gas physical properties (ionization 
fractions and temperature) is computed solving in each cell the 
appropriate discretized differential equations each time the cell is 
crossed by a packet. The reader is referred to Maselli, Ferrara & 
Ciardi (2003) and Maselli, Ciardi & Kanekar (2008) for more de- 
tails. 



2.2 MCLya: line radiative transfer 

MCLya is a numerical scheme for Lya line radiative transfer, 
whose implementation is based on the basic structure of CRASH 
. MCLya in fact uses the same MC sampling and ray-tracing tech- 
niques and it allows for arbitrary 3D hydrogen plus dust density 
distributions, as well as for arbitrary ionization, temperature and 
velocity fields. 

There are three physical processes, included in the code, 
which affect the propagation of the line radiation: Lya line scatter- 
ing, dust absorption and dust scattering. For the sake of simplicity, 
in this paper we concentrate solely on the effect of Lya line scatter- 
ing and we defer the treatment of the interaction between radiation 
and dust to future work. Here we describe the basic structure of the 
algorithm in the absence of dust. For a more complete and accurate 
description the reader is referred to the original paper (VSM06). 

Lya is the strongest HI transition, for which the cross-section 
assumes large values at frequencies near the line center, uq = 
2.466 x 10 lj Hz. It is convenient to introduce the frequency shift: 



v - vq 



(2) 



where Ai/d = (Vth/ 0)1*0 corresponds to the Doppler frequency 
width and Vth is the velocity dispersion of the Maxwellian distribu- 
tion describing the thermal motions, i.e. Vth = (2k bT '/run) 1 ^ 2 = 
12.85 T^ 2 km s _1 , with T4 being the gas temperature in units of 



10 K. The other symbols have the usual meaning. Here we neglect 
turbulent motions, but the option is available for their inclusion. 

The Lya line radiation field is reproduced by emitting photons 
from each source and by following their path through the assigned 
gas distribution until they escape from the simulation box. The lo- 
cation of interaction between the Lya photons and the gas is deter- 
mined by MC sampling the PDF for the line optical depth a photon 
crosses before being scattered, P(~n) = 1 — e~ T ' . In other terms, 
the location of interaction is determined as the cell at which the to- 
tal optical depth from the emission location, t; = ^ . t\ (where 
the sum extends over all the cells crossed by the photon), becomes 
larger than T sca tt ~ — ln(l — £), where f is a random number 
extracted in the interval [0 : 1 [. 

The next step, after assessing the absorption location, is to de- 
termine the photon frequency following a scattering with a hydro- 
gen atom. To do this, the code first converts the frequency of the 
photon from the external (observer) frame, v bs, to the one comov- 
ing with the fluid, v corn , performing a Lorentz transformation: 



'Sobs ^1 



kin " V 



(3) 



where kj„ is the incoming photon direction and V the bulk velocity 
of H atoms. Due to the thermal motion of H atoms, scattering in 
the fluid comoving frame is not perfectly coherent. Within the co- 
moving frameworlQ and neglecting the recoil effect, partially co- 
herent scattering can be described with a simple relation between 
the incoming, Xi„, and the outcoming, x ut, frequency (Dijkstra, 
Haiman & Spaans 2006): 



Va ' kin , V c 



Va 



V th 



(4) 



In the above equation V a is the atom velocity, while ki„ and k out 
are respectively the incoming and outcoming propagation direction. 
The code can model both isotropic and dipolar angular redistribu- 
tion; in this paper we use only the isotropic redistribution and sam- 
ple randomly the outcoming propagation direction. 

Once a new direction and frequency are assigned to the scat- 
tered photon, a new random £ is extracted to determine the next 
scattering location. 

This scheme is repeated until the photon escapes the simula- 
tion volume. 



3 CRASHa: COMBINING CONTINUUM AND LINE 
TRANSFER 

In this Section we describe CRASHa , the first numerical scheme 
which combines the treatment of continuum and line transfer ra- 
diation. As mentioned in the introduction, the algorithm has been 
developed as an extension of CRASH , which provides the treat- 
ment of the ionizing radiation as described in the previous Sec- 
tion and references therein. The extension indeed consists in a 
new algorithm developed to follow the propagation of Lya photons 
through a given gas configuration while it is changed by ionizing 
radiation. In fact, although the continuum photon propagation pro- 
ceeds undisturbed by the Lya radiation field, Lya radiative transfer 
is strongly affected by the change in the ionization state of the gas. 

In order to perform the coupling, it is necessary to intro- 
duce the time evolution for Lya propagation, a feature commonly 

1 Here and in the following we omit the comoving suffix for the sake of 
keeping an easily readable notation. 
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neglected in line radiative transfer codes like MCLya. This is a 
crucial aspect because, due to the resonant scattering nature of 
Lya transfer in a neutral medium, Lya radiation can remain trapped 
for a substantial fraction of the simulation lifetime before being 
able to propagate away from its emission site, while the propaga- 
tion time of the ionization front can be much shorter. Thus, the 
change in the degree of ionization affects the propagation of the 
Lya photons, while the latter induces no back reaction on the gas. 
Note that, although Lya photons, via scattering, can transfer some 
of their energy to the gas and heat it, in typical situations the ef- 
fect is negligible and thus such heating is not generally included 
in Lya radiative transfer codes. We defer the investigation of this 
issue in more detail to future work. 

To correctly model the simultaneous propagation of the two 
radiations a combined approach is needed. This is a challenging 
task because of the very different nature of continuum and line 
transfer, in terms of e.g. their path (straight line versus random 
walk) and time-scales (see discussion above). The above differ- 
ences are reflected also in the numerical implementation of line 
and continuum radiative transfer. For example, while in the case of 
ionizing radiation the time needed for a photon packet to travel a 
given distance does not depend sensibly on the physical properties 
of the gas but only on the physical distance crossed, the propagation 
of Lya photons is very sensitive to the ionization state of the gas 
and extreme configurations can be faced, in which the Lya photons 
scatter for the entire simulation time trapped in few cells without 
exiting the simulation volume. 

As the ionizing radiation scheme has not been modified, it will 
not be discussed further and in the following we will focus on de- 
scribing the details of the line transfer part of the algorithm. 

The Lya radiation is discretized in a large number of photons 
whose emission and propagation is dictated by the time-scale at- 
tached to the ionizing radiation evolution. In this way we are able 
to model the change in the Lya propagation due to the variations in 
the gas ionization state. To correctly model the propagation of Lya 
photons we need to follow every single scattering. As this would re- 
quire a very large computational time, we use a statistical approach 
to the Lya treatment. We have compiled 1085 tables by running 
MCLya, in order to describe the physical characteristics of a pho- 
ton after a scattering depending on the temperature and density of 
the gas and on the incoming photon frequency (see Appendix). The 
following part of this Section is dedicated to a description of the 
various steps of the implementation. 



3.1 Emission of Lya photons 

Every Lya emission is characterized by the generation of iV 7i z 
Lya line photons emitted at the same time, t\ m l . The parameter 
N y j is chosen to optimize the resolution and the code performance. 
The code allows for two different methods for photon emission. In 
the first method the emission is regularly spaced in time as in the 
continuum emission. If, as in Section [2~T1 we define N e m,l as the 
total number of emissions of line photons, in this case: 

A ■ tsim 

'em.l — * X T7 ! (?) 
JV em J 

with i = 0, (N emil - 1). 

An alternative criterion for the emission follows the evolution of 
the ionization structure. In this case the emission time, t l em l , is 
linked to the volume averaged H ionization fraction, xmi.em, and 
Lya photons are emitted at the time t z eml when: 



XHII,em = i X AXHII. (6) 

Axhii = i/N ern , t i is the chosen HII fraction variation in the gas 
and i is an integer that covers values between and N em ,i — 1. 
While in the first formulation a constant Lya emission rate is as- 
sured, in this case the emission rate is higher when the ionization 
state of the gas changes faster. In order to reproduce a constant 
emissivity even in the second formulation, we assign a weight to 
each photon emitted at the i-th step: Wp h = {t l em l —t^A /t sim . 
When a Lya spectrum is built, each photon contributes according 
to its weight. This allows to modulate the emission of Lya photons 
based on the change of the ionization degree (and thus to better 
sample the effect of ionization on Lya scattering) and at the same 
time to have a constant Lya photon rate. 

In the following tests the emission is assumed to be isotropic, 
but it is always possible to account for an arbitrary angular PDF. 

Every emitted Lya photon k (k £ [1, 7V 7i ; x N em j]) is de- 
scribed by its frequency in the comoving frame x in ^ (in this case 
we assume a monochromatic spectrum with Xi„ t k ~ 0, but a dif- 
ferent spectrum can be used), position (which coincides with 
the source location), direction of propagation k in ,k, optical depth 
at which the scattering takes place T sca tt.k (as defined in Sec. l2.2t 
and a characteristic time t c h,h — Kmj that is used to evolve the 
photon along the simulation timeline (see next Section). At any step 
of the simulation the fc-th photon is always described by the quan- 
tities {Xi n , p, ki„, T SC att, t c h), where the index k has been omitted 
for clarity. In the following, we will always omit it. 

3.2 Propagation of Lya photons 

In Section [27X1 we have seen how the physical time of the simula- 
tion is driven by the emission of packets of ionizing radiation dis- 
cretized in time units, dt. We are interested now to link the propa- 
gation of a Lya photon to this timeline. 

Let's assume that an ionizing photon packet has been emit- 
ted at tl m c , that the physical state of the gas has been evolved 
between t\ m c and t J e ^f c , and that a Lya photon is emitted at the 
same time f em l = ii m , c ; then its characteristic time is assigned 
the value t c h = t], m l . The propagation of the Lya photon along 
the direction k is followed between t c h and t J e ^ c , and the line opti- 
cal depth encountered along the path, n, is calculated as described 
in Section |2~2l In each cell crossed by the photon we check if a 
scattering takes place, i.e. if r; becomes larger than T 3ca tt- If there 
is no scattering, we follow the propagation until t{^ c and at this 
point we store the photon's frequency Xi n , the updated position 
p = p + (c dt)]i, and characteristic time t c h = tim.c- Propa- 
gation direction k; n and optical depth for scattering T aca tt remain 
unchanged. These information will be used to follow the photon 
evolution in the next time unit. We define this photon as "active", 
in the sense that it is not trapped by scattering inside a cell but will 
resume its propagation in the next time unit. 

Let's consider now the case in which the photon scatters dur- 
ing the time unit. Unlike MCLya, this code does not follow ev- 
ery scattering inside the cell, but determines the properties of the 
outcoming photon by interpolation of pre-compiled tables (see Ap- 
pendix A). Given the gas temperature, T ce a, the line optical depth, 
T ce u, of the cell where the scattering takes place, and the frequency, 
Xi n , of the incoming photon, a linear interpolation of the tables is 
performed to obtain the distribution of frequencies of the outcom- 
ing photon, x out , and of the time interval that the photon is ex- 
pected to spend inside the cell due to scattering, t sca tt. From these 
distributions the code extracts the values for x out and t 3ca tt that 
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will be assigned to the photon. This approach allows for a tremen- 
dous gain in computational speed by adopting a statistical descrip- 
tion of the scatterings that occur to the photon inside a cell, without 
following each one individually. The characteristic time is updated 
as t ch = tch + t aC att- If t ch > tit^ c the photon is put in a "stand- 
by" mode and its propagation is resumed (with a new T ce u) only 
when the simulation time becomes larger than t c h . 

The procedure described above is repeated for (i) all the 
Lya photons emitted at t r em l , (ii) all the Lya photons "active" at 
t? am , c and (iii) all the Lya photons that exit the "stand-by" mode 
in this time unit. Then, a new ionizing photon packet is emitted at 
tim,c and after it has been evolved up to t{^ c the Lya cycle starts 
again: all the "active" photons are evolved from tf££ c to ^„ 2 C ; 
if there are "stand-by" photons with tim, c < tch < ttrn.c they 
are turned into "active" photons and evolved until t 3 e ^ c ; if new 
Lya photons are emitted in this time unit they as well are evolved 
until t|+%. 

If the ionizing radiation crosses a cell in which a Lya photon 
is trapped by scattering, the change in the physical conditions of the 
cell should be taken into account, as this affect the characteristics 
of the outcoming photon. As an example, let's assume that a Lya 
photon scatters in a cell at t 3 ,o and that the time at which it exits the 
"stand-by" mode is t c h = t s ,2- If the ionizing radiation crosses that 
cell at a time i S) i such that t 3 ,o < t S) i < t s .2, the physical con- 
ditions in the cell change. To take into account the effect on x ut 
and tseatt, we recalculate them from the tables by using the val- 
ues T C ell(t a ,i), T C eii(ts,i) modified by the ionizing radiation and 
Xin(t s ,o) as we do not have any information on the frequency of 
the Lya photon at t s ,i. 



3.3 Spectrum of Lya photons 

When photons exit the simulation box, their frequencies are col- 
lected to calculate the outcoming time integrated spectrum. As 
discussed in Section 13.11 each photon is counted according to 
its weight. To show the probability distribution of the outcoming 
Lya radiation, spectra are normalized to the sum of all weights. The 
profile of the final spectrum strongly depends on the choice of the 
integration time. To build a spectrum we define an initial, t l out , and 
a final, t^J, time. All the photons that escape from the box in the 
interval ]t l out ; t^] will contribute to the spectrum of the source. At 
the next output, all the photons collected in the interval ]t^J ; t^] 
will be used to build the spectrum. This procedure is followed until 
the end of the simulation. As in the case of the emission described 
in Section [3~Tl the spectra can be produced regularly spaced in time 
or linked to the evolution of the ionization structure. Thus, the time 
of the outputs is regulated by equations [5] and [6] where t\ m ; and 
N enl} i are replaced by t l out and N out , respectively. 

Spectra can also be built by choosing a pre-determined line of 
sight. 



4 RESULTS 

In this Section we perform tests for the parallel propagation of 
ionizing and Lya radiation, which show how the evolution of the 
ionization structure alters the Lya spectra of the outcoming radi- 
ation. All the tests have the same initial conditions, unless stated 
otherwise. We use a simulation box of 30 pc on a side, divided 
in 128 3 cells. A monochromatic ionizing source, emitting photons 
with energy equal to 13.6 eV, is located at the center of the box; the 



ionizing photon rate is 5 x 10 49 s _1 . The ionizing radiation is dis- 
cretized in N p = 10 7 photon packets. The same source emits also 
a Lya monochromatic radiation. As we want to construct spectra 
at a fixed distance from the source, we distribute the gas (H only, 
with density tih = 1 cm~ 3 , A*hi ~ 5 x 10 18 cm~ 2 , and tem- 
perature T = 10 4 K) in a sphere of radius r sp h = 15 pc around 
the central source. Outside the sphere the density is set to zero, so 
that no interaction between radiation and gas takes place. The gas 
is initially neutral. Every simulation is carried out for a physical 
time of tsim = 10 yr. In our reference runs we have N ou t = 50 
outputs and N em ,l = 100 emissions of Lya photons, each with 
N lt i — 10 4 photons. Both the emissions and the outputs are dic- 
tated by the evolution of the ionization field. We will discuss the 
effect of a different choice for N em ,u N^ t i and N ou t at the end of 
the Section. In the following we present the results of our simula- 
tions for different choices of the dynamical state of the gas. The 
spectra shown are obtained integrating on all directions in order to 
achieve a better resolution, given the set of chosen parameters. 



4.1 Static sphere 

In this first test, the gas has no bulk velocity with respect to the cen- 
tral source. The top-left panel in Figure[T]shows the spectra emerg- 
ing from this configuration at times corresponding to volume av- 
eraged ionization fractions xhii = 0.3, 0.5, 0.8, 0.9, 0.99, which 
could be regarded as spectra of the source observed at different 
times elapsed since the source switches on. 

The spectra shown here and in the rest of the paper have been 
built as described in Section 3.3., i.e. collecting the Lya photons 
escaping the system when xhii falls in an interval Axhii = 0.02 
centered on the xhii value selected for the output. The curve cor- 
responding to xhii = 0.99 is instead a collection of the escaping 
photons which starts when xhii = 0.98 and ends at t s im- For the 
combination of parameters chosen for the tests, at xhii < 0.3 the 
number of escaped Lya photons is not sufficient to build a spectrum 
and also for xhii = 0.3 the outcoming spectrum is very noisy. As 
expected, the spectra exhibit the two symmetric peaks character- 
istic of this configuration, although a direct, quantitative compari- 
son with previous works (i.e. Dijkstra, Haiman & Spaans 2005) is 
not possible, as none included the effect of ionizing radiation. As 
the ionization increases, the peaks move towards x — because 
Lya photons encounter less and less HI atoms along their path. At 
the same time, the width of the peaks become smaller. In this sce- 
nario we do not see a spectrum peaked at x = because the gas 
never gets completely ionized and, for a static configuration, also a 
little fraction of neutral gas far from the location of emission has a 
non negligible optical depth for photons in the line center. 



4.2 Expanding and collapsing sphere 

In a more interesting case we simulate a homogeneous spherical 
cloud that collapses or expands. In these tests we sample a velocity 
field, in the gas sphere, described by V(r) = V m axr /r sp h + Vo. 

Initially we choose Vmax = ±200 km s _1 and Vo = 
km s _1 . The resulting spectra extracted at the same ionization 
fractions as for the static case are shown in Figure[T]with a positive 
and a negative value for V ma x (top-right and center-left panels re- 
spectively). As expected, the plots show specular Lya spectra, due 
to the opposite direction of the bulk motion. When the gas is ex- 
panding the outcoming radiation is on the red part of the spectrum, 
while it lays on the blue side when we consider negative values for 
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Figure 1. Spectra of Lyct outcoming radiation at times corresponding to ionization fractions Xhii = 0.3, 0.5, 0.8, 0.9, 0.99 (with the exception of the bottom- 
right panel for which xhii = 0.3, 0.4, 0.5, 0.97). The dynamic condition of the gas is the following: static gas (upper-left panel), homogeneous spherical 
cloud expanding and collapsing with velocity increasing with distance from the source (upper-right and center-left panels), homogeneous spherical cloud 
collapsing with velocity decreasing with distance from the source (center-right and bottom-left panels), shell expanding at constant velocity (bottom-right 
nanp.1V 
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the velocity. This happens because the photons are seen Doppler 
shifted according to the velocity of the atoms. This means that, if 
an atom has a positive velocity, in the atom rest frame a blue pho- 
ton becomes a line center photon and is easily blocked by the higher 
optical depth (compared to the optical depth in the wings). Thus, a 
photon can escape only if it is shifted by scattering to the red side of 
the line. Differently from the static case, here it is possible to have 
Lya radiation at the central frequency x = also when ionization 
is not complete. This is a consequence of the fact that, due to the 
Doppler effect, the line center photons are seen in the atom rest 
frame as red (expanding gas) or blue (collapsing gas) photons, i.e. 
in the wing, and thus encounter a lower optical depth. In addition, 
the higher is the absolute value of the velocity the bigger is the 
shift, so when the continuum radiation ionizes the regions closer 
to the source (which have lower velocity and as a consequence a 
higher contribution to the opacity), the optical depth at the center 
decreases significantly because the external neutral layers of gas 
(with a higher velocity) give only a minor contribution. As a result, 
as ionization proceeds, we start seeing an increasing emission at 
the central frequency. More specifically, the spectrum correspond- 
ing to xhii = 0.3 shows radiation at x = and a residual in the red 
(blue) part of the spectrum for positive (negative) velocities. As ion- 
ization proceeds, the residuals become less pronounced and move 
towards the center, while the central radiation becomes stronger. 
In the last spectrum, when the gas is 99% ionized, there is still a 
residual because of the remaining neutral hydrogen fraction in the 
most distant regions of the gas sphere, where the photo-ionization 
rate is suppressed by geometrical dilution and by the residual inner 
opacity. 

In the third case we consider a gas sphere collapsing with in- 
creasing velocity towards the center. The center-right panel of Fig- 
ure Q] shows Lya profiles generated with a bulk motion character- 
ized by V m ax = 200 km s^ 1 and Vo = —200 km s _1 . These 
spectra are very similar to the ones obtained in the previous case, 
but as the absolute value of the velocities increases towards the cen- 
ter, the residual blue part of the spectra is more spread. Note that 
also in this case the residual is reduced as ionization proceeds and 
the contribution to the opacity from the inner layers of gas is sup- 
pressed. 

In the last case (Fig. [T] bottom-left) we consider a gas which 
is collapsing near the source while the outer shells are expanding, 
with V m ax = 200 km s" 1 and V = -100 km s -1 . The first 
Lya spectrum (at xhii = 0.3) is dominated by photons in the 
blue part and just a residual is present on the red side. In fact, blue 
photons have a larger probability to escape because the ionization 
front has not yet propagated far enough to suppress the contribution 
to the Lya gas opacity from the gas collapsing towards the source. 
As the front proceeds ionizing the neutral hydrogen with negative 
velocities (spectra at increasing Xhii), the red part of the spectrum 
becomes stronger while the blue part is suppressed, until, in the 
configuration with Xhii = 0.99, it is smaller than the red one. As 
in the other tests, the increment in the ionization degree reduces 
the number of scatterings, with the consequence of moving the two 
peaks towards the center, increasing their height and reducing their 
width. 

4.3 Expanding shell 

In this Section we examine the case of a Lya source surrounded 
by an expanding shell, which has been extensively studied also 
by other authors (Ahn, Lee & Lee 2004; Hansen & Oh 2006; 
VSM06). Here, we are interested in the effects introduced on the 



Lya spectrum by an ionizing source inducing a time evolution 
of the neutral gas in the shell. To simulate this configuration we 
have chosen a homogeneous density nu = 15 cm -3 , temperature 
T = 10 4 K, and radial velocity V = 300 km s _1 . All the gas is 
distributed within a shell of thickness 4 pc located at a distance of 
10 pc from the source, while no gas is present outside the shell. 
The corresponding column density is Nm ~ 2 x 10 20 cm -2 . To 
show the Lya spectrum time evolution, we choose to plot the pro- 
files corresponding to ionization fractions in the shell of xhii = 
0.30, 0.40, 0.50, 0.97 (Fig.Q] bottom-right panel). For a better un- 
derstanding of the spectral features, it is useful to discuss the pos- 
sible different paths for an outcoming photon. Following VSM06 
(see their Fig. 12), we divide the outcoming photons in three differ- 
ent groups, depending on their scattering history: 

• Backscattered photons: photons that, after scattering in the 
shell, travel inward across the empty space before crossing again 
the shell. As these photons undergo multiple scatterings with the 
gas, they can escape once they are shifted on the red side of the line 
where the optical depth of the expanding shell is smaller. 

• Diffused photons: all the photons which are diffused in the 
shell until they escape without backscattering. We expect these 
photons to contribute to a red bump in the spectra whose shift from 
the line center and intensity will depend on the neutral gas den- 
sity and on the shell velocity. Typically the frequency shift will be 
smaller than for backscattered photons as the number of scatters 
before escape is on average lower. 

• Directly escaped photons: photons that have no interaction 
with the gas and keep their initial frequency; in our case this group 
of photons will produce a peak at x = 0. 

It is important to underline that every group has a different char- 
acteristic time for escaping. In fact, directly escaped photons travel 
the shortest path. On the other hand, diffused photons scatter in 
a volume smaller than the backscattered photons and thus escape 
faster; therefore the red bump associated to the diffused photons 
will typically appear before the feature produced by the backscat- 
tered photons. 

Keeping in mind all the possible paths for Lya photons, let us 
analyze the features in the spectra shown in the Figure Q] (bottom- 
right panel). The first profile (xhii ~ 0.3) exhibits a peak on the 
red side of the central emission, due to the diffused photons, that, 
as already mentioned, escape faster than backscattered photons and 
can already been seen in the initial stages of the shell ionization. 
Directly escaped photons are present as well and their abundance 
increases with time. In the profile corresponding to an ionization 
degree of xhii = 0.4 we can clearly see that a secondary bump 
is forming at lower frequencies. At this stage of the evolution, the 
backscattering photons are starting to escape from the shell with 
a frequency that is more shifted respect to the other photons, as 
explained above. When the ionization degree is xhii = 0.5, the 
backscattering bump is visible and dominant on the red peak due to 
diffused photons. In the profile corresponding to xhii ~ 0.97 the 
ionization front has suppressed most of the neutral gas and only a 
negligible fraction of Lya radiation interacts with the residual gas 
in the shell; the result is a small fraction of photons shifted on the 
spectrum's red side. 

4.4 Effect of ionizing radiation 

In the previous tests we have discussed how our time dependent 
treatment of the Lya radiation allows to correctly establish the ap- 
pearance at different times of spectral features which are usually 
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integrated in the emergent spectra predicted with time independent 
formulations. Here we investigate further on the importance of the 
joint propagation of Lya and continuum radiation, by comparing 
results from two different approaches to simulate the Lya spectra. 
The first one, widely used in the literature, performs a Lya radiative 
transfer on a gas configuration given as initial condition, which can 
be e.g. a constant density field or a snapshot of a numerical simula- 
tion. In this case, the gas configuration is kept constant throughout 
the entire Lya radiative transfer and the Lya spectra are built once 
all the Lya photons have escaped the simulation volume. The other 
approach is the one described in this paper, i.e. starting from an 
initial gas configuration, the parallel propagation of continuum and 
line photons is followed and the Lya spectra can be built at different 
times taking into full account the changes in the Lya propagation 
due to the variations in the neutral gas distribution. 

To show the impact of the two different approaches on the out- 
coming Lya spectra, we consider the same configuration described 
in Section |4~3l i.e. an initial neutral expanding shell which is ion- 
ized by a central source emitting also Lya photons. We compare the 
spectra obtained with CRASHaat the times when the gas configu- 
ration is characterized by a mean ionization fraction inside the shell 
of xhii = 0.32 and 0.66, to those obtained running MCLya on the 
same gas configurations. While with MCLya the spectra are built 
by integrating over the Lya photons once they have all escaped the 
fixed HI distribution, the CRAS Ha algorithm allows to account for 
the impact on the emergent spectra of the ionization history which 
led to those configurations. The results are shown in Figure|2] 
A substantial difference is clearly visible in the spectra correspond- 
ing to xhii = 0.32. The Lya spectrum simulated by MCLya 
is characterized by two bumps associated with the diffused and 
backscattered photons. A very different profile is obtained with 
CRASHa , where no backscattered photon has escaped at this time 
and only the single red peak corresponding to the diffused pho- 
tons is present. As already underlined in Section l4~3l the absence of 
backscattered photons is due to their larger escaping time. The pro- 
files corresponding to xhii = 0.66 are much more similar, because 
at this stage a significant fraction of the backscattered photons had 
enough time to escape the shell. Nevertheless, there is still a differ- 
ence in the amplitude of the peak at x — and of the bumps from 
the backscattered photons, which are also slightly more shifted. 
This difference is due to the memory of Lya photons emitted in 
the previous stages of the source activity, when xhii < 0.6. In the 
MCLya treatment all the emitted Lya photons see the same mean 
shell opacity and the probability to have a direct escape is signif- 
icantly higher then in the CRASHa run, in which the Lya photons 
see on average a larger shell opacity. As a consequence, the fraction 
of photons with x = in CRASHa is smaller and, at the same time, 
the photons that have remained trapped for a longer time exhibit 
a larger amplitude of the spectra and a larger shift to the red side 
of the central frequency. The time elapsed between the two spectra 
considered above is only about 100 yr. The time interval in which 
deviations from the "instantaneous picture" are significant is thus 
too small to allow observations to capture these stages. However, 
the example is useful to illustrate the relevant features and advan- 
tages of our approach. Furthermore the deviations found could af- 
fect the gas state, e.g. its spin temperature, independently from their 
observational detectability in the spectra. 

A similar test has been performed to infer the observability of 
the deviations in the spectra on larger scales. In this case we use a 
simulation box of 200 kpc on a side and a photo-ionization rate of 
the central source N 1 = 10 s . In this test the gas is distributed 
in a sphere of radius r sp h = 70 kpc, with density nu = 0.01 cm -3 




x 



Figure 2. Comparison of Lya spectrum obtained with MCLya and 
CRASHa applied to the same gas distribution of an expanding shell. For 
the application of MCLya a fixed gas configuration is used, while in our 
approach all the evolutionary stages are taken into account (see text for de- 
tails). 



(corresponding to Nhi ~2x 10 21 cm 2 ), while we keep the tem- 
perature T — 10 4 K; we also assume a velocity field corresponding 
to a Hubble expansion with H —190 km s _1 Mpc -1 . The simula- 
tion is carried out for a physical time of t 3 i m = 10 8 yr. In this case 
(Fig. |3j the first spectrum is captured at the time t = 4.0 x 10 6 yr 
(corresponding to a mean ionization fraction of xhii = 0.24) and 
the second at the time t = 8.3 x 10 B yr (xhii = 0.40). While the 
MCLya profile in the left panel of Figure [3] is characterized by a 
large red bump and a small blue one, with the CRASHa approach 
we find a lower fraction of photons escaping with blue frequencies 
and a larger red bump which is slightly shifted on redder frequen- 
cies. This effect results from keeping memory of the ionization his- 
tory, since we take into account that before reaching the observed 
gas configuration most of the Lya photons have been trapped on 
the boundary of the growing ionized region. The integrated optical 
depth along the full path traveled before escaping is therefore larger 
respect to the one computed in the MCLya approach. This effect is 
also visible in the right panel of Figure [3] where a larger red shift 
is present. In this case a larger fraction of trapped photons are now 
free to escape and the shift is more evident. 

The major result of these tests is that the emerging spectra 
keep memory of the ionization history which generates a given ob- 
served configuration and, to properly account for this effect, the 
self-consistent joint evolution of line and ionizing continuum radi- 
ation followed by our scheme is necessary. The extent of the dif- 
ference between the two methods depends on the particular case 
considered, but it can be substantial and can thus affect the physical 
interpretation of the problems at hand. In a forthcoming study we 
will investigate in more detail which are the objects/configurations 
for which the ionization effects are expected to be relevant in shap- 
ing the observed Lya spectrum. 

In addition, the time evolution that builds up the Lya radiation 
field can be important when calculating the impact of such radiation 
on gas properties like the spin temperature, which is relevant for 
the prediction of the observability of 21 cm emission from neutral 
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hydrogen at high redshift. We plan to include the computation of 
these effects in a forthcoming extension of the code. 



4.5 Dependence on input parameters 

In this Section we investigate the dependence of the final spectra on 
the parameters N em ^, N Jt i and N out defined in Sections |3, II and 
13.31 The convergence tests presented are performed by adopting the 
same conditions of the static case (see Sec. 4.1). 

First we compare spectra built by integrating over the full sim- 
ulation time, i.e. we set N out = 1. We have done several runs vary- 
ing the number of Lya emissions (iV em ,; — 2, 10, 50, 100) and 
setting Njj = 10 4 . Lya photons are emitted according to equa- 
tion[6] i.e. at regular ionization intervals. For N em ,i = 2, the two 
emissions are performed at the beginning and once the ionization 
degree has become stationary. The results are shown in the top-left 
panel of Figure [4] From an inspection of the Figure it is clear that 
2 emissions are not accurate enough to properly describe the emer- 
gent Lya spectrum. This is a consequence of the method described 
in Section 3.1 for weighting the contribution to the spectrum from 
photons emitted at different steps. In fact, as the time between the 
two emissions is large, the weight assigned to the photons of the 
second (and last) emission is so large that the contribution from the 
photons of the first emission is negligible. Moreover, at the time of 
the second emission ionization is almost complete and as a conse- 
quence the spectrum has two thin peaks very close to the central 
frequency and shows no large frequency shift due to scattering in a 
gas with larger opacity. So, in this case the outcoming spectrum is 
the result of Lya transfer in an almost ionized sphere and the stages 
through which the gas reached such configuration are completely 
neglected. As the number of emissions is increased, the accuracy 
improves and convergence is reached when N em ,i > 50. 

We have then checked convergence of instantaneous (vs in- 
tegrated) spectra, setting N out = 10 and performing again four 
runs with different numbers of Lya emissions (N em .i = 2, 10, 50 
and 100) and iV 7l ; = 10 4 . Figure [4] displays two sets of spectra, 
each corresponding to a fixed xhii value: 0.4 (top-right panel) and 
0.7 (bottom-left panel). Effectively, the spectrum corresponding to 



N ern ,i = 2 is built only with photons from the first emission, as the 
second is performed when xhii > 0.7. At xhii = 0.4 the spec- 
tra are all very similar because the gas opacity is sufficiently high 
to trap Lya photons and those emitted close to the time when the 
spectrum is built are still scattering in the gas. On the contrary, as 
ionization proceeds, photons emitted at later times encounter less 
neutral hydrogen and escape more easily, and their contribution to 
the outcoming spectrum increases. As a consequence, the case with 
Nem,l = 2, when Lya photons have been emitted only at the begin- 
ning of the simulation, displays a broader and more shifted profile. 

Finally, we use different values of /AT 7|1 = 10 3 , 10 4 , 10 5 , with 
N em j — 100 and N out = 50, to determine the smallest number 
of Lya photons needed in each emission in order to achieve conver- 
gence. In the bottom-right plot we show Lya profiles corresponding 
to an ionization fraction of Xmi = 0.7. As expected, the increas- 
ing number of Lya photons in every emission produces a smoother 
profile. We find that, for the configuration considered, the emis- 
sions characterized by 10 3 photons are not accurate enough, and at 
least 10 4 photons per emission are needed. 



5 SUMMARY 

In this paper we have presented CRASHa , the first radiative transfer 
code for cosmological application that follows the parallel propa- 
gation of Lya and continuum photons. Since Lya propagation is 
dominated by resonant scattering with neutral gas, the effect of a 
rapid change in the degree of gas ionization can affect the features 
of the emerging Lya spectrum. To investigate this issue, we have 
developed in the continuum radiative transfer code CRASH (Cia- 
rdi et al. 2001; Maselli, Ferrara & Ciardi 2003; Maselli & Ferrara 
2005, Maselli, Ciardi & Kanekar 2008) a new algorithm to follow 
the propagation of Lya photons through a gas configuration while 
it is changed by ionizing radiation. 

In order to perform the implementation, it has been necessary 
to introduce the time evolution for Lya propagation, a feature com- 
monly neglected in line radiative transfer codes. This is a crucial as- 
pect because, due to the resonant scattering nature of Lya transfer 
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in a neutral medium, Lya radiation can remain trapped for a sub- 
stantial fraction of the simulation lifetime before being able to 
propagate away from its emission site, while the propagation time 
of the ionization front can be much shorter. Another challenge of 
the implementation has been to reduce the computation time for 
the Lya scattering. In fact, to correctly model the Lya propagation 
every single scattering should be followed. As this would require 
prohibitively large computational times, we have used a statisti- 
cal approach to the Lya treatment. We have compiled tables us- 
ing MCLya (Verhamme, Schaerer & Maselli 2006) to describe the 
physical characteristics of a photon escaping from a gas cell where 
it was trapped by scattering as a function of the temperature and 
density of the gas as well as of the incoming photon frequency. The 
tables are called within CRASHa . With this statistical approach we 
experience a drastic reduction of the computational time and, at 
the same time, an excellent agreement with the full Lya radiative 
transfer computations. 



We have discussed the details of the code implementation and 
tested it for several gas configurations, including static spherical 
gas distribution, expanding and collapsing spheres and expand- 
ing shell. For all the configurations analyzed, CRASHa reproduces 
emerging spectra with the qualitative features expected from the- 
oretical models and discussed previously in the literature, while a 
more quantitative comparison has not been feasible as CRASHa is 
the first code which couples the continuum and line transfer. With 
this implementation, it has also been possible to investigate how 
the line shape of the emergent spectra evolves with the gas ioniza- 
tion. Although the specific results depend on the geometry of the 
gas and on the velocity field, common trends are found. The main 
results can be summarized as follows. 

• While ionization proceeds the peaks on the blue/red side of 
the line center move closer to the central frequency, getting thinner 
and higher, as expected for a gas configuration with progressively 
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decreasing optical depth. Depending on the gas configuration (e.g. 
in case of an expanding shell), more complex features arise that can 
be associated to the different paths followed by the Lya photons 
before escaping. 

• The emerging spectra keep memory of the ionization history 
which generates a given gas configuration. 

• The novel approach to Lya transfer developed in 
CRASHa allows to resolve the emergence of different spec- 
tral features at different times during the evolution of the ionization 
field. Features emerging on different time scales are typically 
associated to the various paths traveled by the photons before 
escaping. 

• A comparison between our new algorithm to follow the prop- 
agation of Lya photons and a full line radiative transfer shows an 
excellent agreement for different gas configurations and an enor- 
mous gain in computational speed. 

In order to account for the effects discussed above a self- 
consistent joint evolution of line and ionizing continuum radiation 
as implemented in CRASHa is necessary. A comparison between 
the results from our code and from Lya scattering alone on a fixed 
density field shows that the extent of the difference between the 
emerging spectra depends on the particular configuration consid- 
ered, but it can be substantial and can thus affect the physical inter- 
pretation of the problem at hand. 

A detailed discussion on which are the objects/configurations 
for which the ionization effects are expected to be relevant in shap- 
ing the observed Lya spectrum is deferred to a forthcoming publi- 
cation. Nevertheless, we have here discussed two specific test con- 
figurations for which the coupling of continuum and line radiation 
would be necessary in order to recover correctly the emergent pro- 
file. These differences are due to the time evolution feature intro- 
duced in CRASHa for Lya photons which allows to keep track of 
the ionization history imprint on Lya profiles for the first time in 
the literature. 

The time evolution that builds up the Lya radiation field can 
be furthermore important when calculating the impact of such radi- 
ation on gas properties like the spin temperature, relevant to predict 
the observability of 21 cm radiation from the early universe. In a 
forthcoming extension of the code CRASHa , we plan to include the 
self-consistent calculation of the impact of the Lya radiation on the 
gas temperature, together with the contribution to the Lya radiation 
from recombinations occurring in the gas. 
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APPENDIX A: TABLES FOR LYa SCATTERING 

Following the single scatterings of each Lya photon can be ex- 
tremely expensive; in order to avoid it we have built tables by using 
a statistical approach that allows to retrieve the frequency of the 
outcoming photon, x ou t, and the time interval for which the photon 
is trapped in the gas by the scatterings, t 3ca tt, given the frequency 
of the incoming photon, Xi n , the gas temperature, T ce u, and optical 
depth, T ce u, of the cell where the scattering takes place. We adopt 
the opacity at line center to characterize the optical depth of the 
gas in a cell. Notice that the frequencies are always meant in the 
comoving frame. 

The tables are compiled in the following way. Given a value 
for the input parameters (T ce u, T ce u and Xi n ), we run the MCLya 
code several times (the results converge when 10000 Lya photons 
are emitted) to obtain values for x out and t sca tt which are then 



12 M. Pierleoni, A. Maselli and B. Ciardi 




,10 



X 



Figure Al. Distribution of photon outcoming frequencies generated by Figure A2. Distribution of the time that a photon with Xi n = spends in a 

multiple scatterings inside a cubic cell with T ce u = 10 4 K and r ce ; j = cubic cell with T ce ;; = 10 K and an optical depth r ce ;; = 10 4 (solid line), 

10 6 for an incoming photon frequency of ij„ = 20 (solid line), 30 (dashed 10 5 (dashed line) and 10 6 (dashed-dotted line), 
line) and 40 (dashed-dotted line). 



binned in distribution functions. The x out distribution is binned 
using regular intervals of 0.5. For tscatt we adopt equally spaced 
logarithmic bins of 0.014. It is important to note that the value of 
tscatt, which depends on the distance traveled by the photon, is 
linked to the size of the cell. Thus, the tables are compiled for a ref- 
erence cell size, d c , r ef, but anytime they are accessed by CRASHa , 
the value of t scat t obtained needs to be rescaled for the actual cell 
dimension, d c . The ranges that the tables cover are: 

• Temperature: 10 K < T < 10 5 K 

• Optical depth: 1 sC r < 10 6 

• Frequency: —100 < x < 100 

An example of x out and tscatt distributions is given in Fig- 
ures I Al I and IA2I Figure IA1I shows the distribution of x out for 
T ce u — 10 4 K, T ce u = 10 6 and different values of Xi„ (xi n =20, 
30 and 40). Figure lA2l shows how the t scat t distribution changes for 
different values of the optical depth, r ce u = 10 4 , 10 5 , 10 6 ; here we 
fixed T ce u — 10 K and x in = 0. 

The above tables are accessed by CRASHain the following 
way. If T ce n, T ce u and Xi„ are within the range covered by the 
tables, the distributions for tscatt and X ou t <H6 calculated by a lin- 
ear interpolation and then the value used in CRASHafor tscatt and 
Xout is obtained by MC sampling the interpolated distributions. 
The interpolation scheme for tscatt is as follows. As a first step, 
the values of temperature Ti and T2 closest to T ce ii for which 
T\ ^ T ce u < T2 are found. The same is done for the optical 
depths n and T2, and the frequencies xi and x%. The weights to 
be assigned to each value are derived by linear interpolation. As an 
example we can consider a simple ID case, with linear interpola- 
tion only on temperatures. In this case, t sca tt = wri tscatt (Ti) + 
Wt 2 t S catt (T2), where the weights are w Tl = (T 2 - T ceH )/AT 
and wt 2 = {T cM - Ti)/AT, with AT = T 2 - Ti. The same 
procedure, extrapolated in 3D, produces a distribution of t sca tt that 
will be randomly sampled in CRASHa . The interpolation for x ut 



follows the same steps, after a shift in the frequency space has been 
performed to assure that the correct distribution is obtained and 
no spurious peak forms. To understand how the x out interpolation 
works an easy example can be useful. Let us assume that our pur- 
pose is to reproduce the dashed profile centered in x = 30 (Fig. lAU 
starting from the two profiles at x = 20 and 40. If the interpolation 
were performed bin by bin without any previous shift (e.g. the bin 
[19.5-20[ of the solid curve with the bin [39.5-40[ of the dashed- 
dotted curve, and similarly for all the bins), the result would be two 
smaller peaks centered at x = 20 and 40, rather than one centered 
at x = 30. To perform a correct interpolation we need to center the 
solid and dashed-dotted profiles on x — 30 and then interpolate. 

To check the validity of our approach, we have compared the 
results from the full radiative transfer treatment (using MCLya) 
with a case in which the tables were used. In Figure |A3l the distribu- 
tion in frequency of 10000 Lya photons escaping from a gas with 
temperature T = 8000 K, optical depth r = 10 6 and frequency 
x — is shown. The dotted (solid) line in the upper panel indicates 
the results for the approximate (full radiative transfer) treatment us- 
ing tables built with 10000 photons. In the bottom panel the differ- 
ence between the two distributions is plotted, showing an excellent 
agreement, which is found also for different initial conditions. We 
perform the same check with a more complex gas distribution using 
the expanding shell described in Section |4~3l As the final spectrum 
is expected to have a wider frequency distribution, this test is per- 
formed to check the accuracy of the x out interpolation at larger 
frequency shifts. The results are shown in Figure lA4l for both pro- 
files we have used 400000 photons. Also in this case an excellent 
agreement is found. 

As the changes in the gas properties during a simulation can 
be drastic, sometimes the values r ce ;; and Xi„ can fall outside the 
range covered by the tables (we do not expect T ce u to fall out- 
side the range). In these cases we perform an extrapolation of the 
existing tables. More in particular, for \\xi n \\ > 100 we use the 
same distributions derived for \\xi n \\ = 100. This is a good ap- 
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The use of the tables allows an enormous gain in computa- 
tional speed. 



Figure A3. The upper panel shows the predicted Lya frequency distribu- 
tion for monochromatic line radiation escaping from a gas with temperature 
T = 8000 K and optical depth r = 10 6 . The dotted (solid) line indicates 
the results for the approximate (full radiative transfer) treatment using ta- 
bles build with 10000 photons. The difference between the two distributions 
is shown in the bottom panel. 
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Figure A4. As Fig. lA3l but for a gas configuration resembling a shell with 
temperature T = 10 4 K and column density N m = 2 x 10 20 cm -2 
expanding with an uniform radial velocity V = 300 km s — x . The dotted 
(solid) line in the upper panel indicate the results for the approximate (full 
radiative transfer) treatment using tables build with 400000 photons. 



proximation as at these frequencies the cross section is small and 
Lya scattering rare. The extrapolation for the optical depth works 
differently. If T ce u < 1 no interaction takes place and the photon 
propagates freely. If r ce n > 10 6 we divide the cell into 2 3n sub- 
cells (where n is the number of divisions performed) until each 
sub-cell has an optical depth < 10 B . At this point, every sub-cell is 
treated as a single cell. 



